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We investigate the effect of microbubbles on Taylor-Couette flow by means of direct 
numerical simulations. We employ an Eulcrian-Lagrangian approach with a gas-fluid 
coupling based on the point-force approximation. Added mass, drag, lift, and gravity are 
taken into account in the modeling of the motion of the individual bubble. We find that 
very dilute suspensions of small non-deformable bubbles (volume void fraction below 
1%, zero Weber number and bubble Reynolds number < 10) induce a robust statistically 
steady drag reduction (up to 20%) in the wavy vortex flow regime (Re = 600 - 2500). The 
Reynolds number dependence of the normalized torque (the so-called Torque Reduction 
Ratio (TRR) which corresponds to the drag reduction) is consistent with a recent series 
of experimental measurements performed by Murai et al. (J. Phys. Conf. Ser. 14, 143 
(2005)). Our analysis suggests that the physical mechanism for the torque reduction in 
this regime is due to the local axial forcing, induced by rising bubbles, that is able to 
break the highly dissipative Taylor wavy vortices in the system. We finally show that the 
lift force acting on the bubble is crucial in this process. When neglecting it, the bubbles 
preferentially accumulate near the inner cylinder and the bulk flow is less efficiently 
modified. 



1. Introduction 

Drag reduction induced by injection of small concentration of gas bubbles in a liquid 
flow has b een addressed in sever al physical systems, from turbu le nt boundary layer on a 
flat plate dMadavan et all (1 19841) . iFerrante k, Elghobashl (|2004l ). ISanders etall (|2006h ). 



to channel flows ( Xu et all (|2002h . iLu et all (|2005l )). Despite the efforts dedicated to the 



subject, the problem has not been fully clarified from a fundamenta l physical point of 



view. Different mechanisms suc h as effect i ve com pressibility of the flo w (IFerrante &: Elghobashi 
( 2004 )). bubble deformabil i tv (ILu et all (|2005f) ). or compressibility (|Lo et aLl (|2006f )) . or 



splitting ( Meng &: Uhlmanl ( 1998f ). have been proposed as relevant. A misleading aspect of 



the drag-reduction problem arises from the fact that, depending on the system considered, 
either transient, or spatially dependent (as a function of the bubbles' injection point), 
or statistically steady effects can be observed. To overcome such a diffic ulty, more re- 
cently experiments have been conducted in a Taylor-Couette (TC) set-up, ( Dieridi et all 
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( 20041 ) . Ivan den Berg et~al\ (|2005h . lMurai et all (|2005l ). Ivan den Berg efaH (|2007t )). The 



TC system has two advantages. First, since it is a closed system, exact global energy 
balance relation can be explicitly obtained and drag variations can be evaluated in terms 
of a single globally averaged quantity, the torque. Second, statistically stationary states 
can be reached easily. In order to quantify the level of drag reduction, it is convenient to 
consider the so called Torque Reduction Ratio (TRR) coefficient, viz. 

TRR = 1 — ^ , (1.1) 

where Tb stands for the measured torque in the two-phase system, i.e., with the bubbles 
included and T for the torque in the single-phase flow at the same Reynolds number. 
In the highly turbulent TC system (Re > 10 5 ) drag re d uctio n is mainly attrib uted to 
the bubble deformation mechanism (Ivan den Berg et~al\ (|2005l ). lLu eTaU (|2005f) ) in the 



boundary layer (|van den Berg et all (|2007l )). 

Also at moderate flow-Reynolds numbers i n the range Re = 6 00 - 4500 (the so called 
wavy/modulate-wavy Taylor vortex regime) iMurai et al. (|2005h reported that bubble- 



induced drag reduction up to 25% can be observed for tiny bubble concentrations, 
O(0.1%) in volume. The experimental apparatus considered by these authors is forced by 
rotation of the inner cylinder, while the outer one is kept fixed. The cylindrical enclosure 
is filled by a Silicon oil roughly five times more viscous than water (y = 5- 10 -6 m 2 /s) and 
air bubbles are injected into the fluid. Their typical radius is a = 250 - 300 /im and We- 
ber number We = 2apv^(r~ 1 < 0.6 (where p and a are respectively the density and the 
surface tension of the oil and vt the bubble terminal velocity in still fluid). The bubble 
Reynolds number, based on vt, is Reb = 2a vt/v ~ 7. Therefore bubbles can be regarded 
with good approximation as mono-disperse non-deformable spheres, whose scale is of the 
same order as the smallest scales of fluid fluctuations. The experimental results of Murai 
et al. showed that the drag reduction decreases with increasing Reynolds number. TRR 
was almost vanishing at Re ~ 4000 and became negative (drag increase) for even larger 
Reynolds numbers. Murai et al. explained this drag reduction effect through the bubble 
interaction with the coherent Taylor vortices, producing a vertical elongation of their 
arrangement. 

The features of the Murai et al. experiment, in particular the low large-scale Re num- 
bers and the small (i.e. limited Reb) non-deformable bubbles employed there, make this 
experiment accessible to a comparison with numerical simulations, which we will perform 
in this paper. 

The aim of the present work is to shed more light on the physical mechanism of drag 
reduction observed in the flow conditions of Murai et al. experiment. However, given the 
complexity of this two-phase flow system, appropriate modeling of the dynamics is un- 
avoidable and approximations must be introduced. First, we use an Eulerian-Lagrangian 
algorithm. The forces on each point-bubble are modeled through effective drag and lift 
forces; they act back on the fluid through two-way coupling. Second, in order to make 
the simulation computationally feasible, we have to restrict the flow to a smaller domain 
as compared to the experiment and to adopt different, i.e. periodic, boundary condi- 
tions. Third, in order to achieve good convergence, the bubbles are initially positioned at 
random positions throughout the flow. This is obviously different from the experimental 
procedure. In the experiment bubbles are injected from the bottom. The injection sites 
are close to the inner cylinder walls. The bubbles emerge to the free surface at the top 
after traveling across the cylindrical enclosure, whose aspect ratio is much larger than 
in the simulation (cf. the T values in Table [lj. It is finally on the top TC section that 
the radial mean void fraction is measured. The experimenters assumed that the bubble 
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Figure 1 . Sketch of the Taylor-Couette system 

distribution rapidly reaches a fully developed state that does not keep memory of the 
inhomogeneous injection sites. This assumption justifies the different bubble injection 
procedure adopted in our numerics. 

In spite of these unavoidable differences, we think - and we will further motivate in 
the paper - that the main flow features and the resulting drag reduction phenomena 
are reasonably captured by our numerical approach. In particular, it is of relevance 
to note that drag reduction mechanism observed in the flow conditions examined here 
- which we will call pseudo-turbulent drag reduction - turns out to be substantially 
dif ferent from the mechanism that acts under the highly-turbulent conditions explored 
bv Ivan den Berg et all {2005, 20o 3). Under those conditions the bubble deformation had 
turned out to be crucial, see also lLu et all (|2005h . 

The paper is organized as follows: in section[2j we introduce the geometry and the two- 
phase model system used for the numerical study and present a test in order to validate 
the code. In section [3] our numerical results on torque modulation in the two-phase flow 
and bubble distributions are described and discussed in detail. Section [4] contains com- 
ments and conclusions. In the appendix we derive relations connecting globally averaged 
quantities in the TC system. 



2. Taylor-Couette system and numerical scheme 

2.1. Definitions 

In the Taylor-Couette system the fluid is enclosed between two coaxially rotating cylin- 
ders with radius R and Ri, respectively for the outer and the inner one. The height 
of the cylindrical enclosures, denoted with L, is in the vertical direction. We limit our 
investigation to the case where only the inner cylinder is rotating at a constant angular 
velocity £1 while the outer one is at rest, see Fig. [T] Furthermore, no-slip conditions are 
assumed at the walls. It is convenient to adopt the Reynolds number based on the inner 
cylinder velocity U — Ri^l and the gap width d = (R a — Ri) as a control parameter, that 
is 

Re ^ (R -R i )R i n (21) 
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The response of the system is the time averaged torque on the inner cylinder, which is 
expressed as 

T = 2kR 2 1 Lt wi , (2.2) 

where r W i stands for the mean inner-wall shear stress. In cylindrical coordinates (z,r,9) 
the wall shear stress can be computed as 



-pv \ -^ue\ r =R t 



n 



(2.3) 



where ug is the azimuthal velocity component and the over-bar denotes averaging in 
space over concentric cylindrical surfaces and time. Note that an analogous expression 
for the torque can be written in terms of the outer cylinder, since conservation of angular 
momentum implies {■^ue\ r =R )R 2 , — {-^ug\r=B.i — H)iZ|. At increasing angular velocity 
the flow in between the two cylinders experiences several hydrodynamic instabilities. Af- 
ter the dcstabilization of the la minar flow, the formation of horizontal couples of circular 
rolls, known as Taylor vortices (Taylor ( 19231 ) ). is observed. Then, at further increasing 
f2, axial symmetry is lost and one encounters first steady wavy, then time modulated 



rolls, and finally a fully developed turbulent regime is established, see lAndereck et 
(|l986h . 



The relation linking the torque to the Reynolds number and the gap parameter r\ = 
Ri/R cannot be directly computed from the equation of motion, except for the laminar 
(Couette) regime. To de duce under gen e ral conditions the functional relat i on T(Re, 77) is 
still an open problem see Lathrop et al. ( 1992 ). Esser fc Grossmann ( 19961 ). Eckhardt et a 
(|2000t ). iLim fc Tanl (|2004[) . and Eckhardt et all (|2007i ). However, the torque can be re- 
lated to the mean (time and volume averaged (. . .}) energy dissipation rate (s), namely 


Til 



{e) = 



(2.4) 



Thus (e) and T are directly proportional and a decrease in the torque corresponds to a 
reduction of the internal energy dissipation rate in the system. Therefore, the total drag 
reduction can be deduced from the variations of a single global observable as the torque. 



2.2. The model for the two-phase flow 

The direct numerical simulation (DNS) approac h we adopt is ba s ed on the Eulerian- 
Lagrangian algorithm, here in the form used by iMazzitelli etaH (|2003olf6h . The fluid 
phase is described by the Navier-Stokes equation for an incompressible flow (V • u = 0). 
Bubbles counter-act on the flow by a two-way point-like coupling: 



d t u+(u- V)u 



— - + v V 2 u 



(2.5) 



16 = 



47ra3 xt 



(d t u + (u • V) u - g) , 



(2.6) 



where y/n denotes the instantaneous position of the i-th bubble and g the gravitational 
acceleration. The bubbles are tracked as Lagrangian particles with a model equation for 
wh ich added mass, d rag, l ift , and g r avity a r e taken into a c count . Their dynamics is given 
by iMaxev fc Riievl (|l983l ). lAutonl (|l987[) . lAuton et al] (|l988f ). ICliment fc Magnaudet 



f Note that in lLathrop et al] l| 19921 ) and in Ivan den Berg et al\ (|2005l ) this relation contains 
a typo, namely, a wrong pre-factor 2. See also the appendix. 
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Figure 2. Single-phase flow code validation, (lef t) The c ritical Reynolds number Re 



function of the gap ratio n, a n d comparison with iTavioi (ll923|) narrow-gap limi t (n 
iDominguez-Lerma et al\ (| 1 9841 ^1 . lEsser fc Grossman iil i|1999 ) and lDutcher fc Mullerl (|2007n pre- 



as a 
1), 



dictions, (right) Dimensionless torque G = T/(pu 2 L) versus Re for three different values of the 
gap ratio r;, a nd comparison with Wendt's empirical relation Gw = 1.45(i?e 77) 3//2 /(l — f?) 7//4 , 
IWendtl l|1933f ) . The inset shows a compensated plot G/Gw vs. Re. 



(|l999t ). lMazzitelli et all ( j20036T ): 



dy 
dt 
dy 



v. 



(2.7) 



■37 = 3 (d t u - 



V)u)--(v 



u) - 2g - (v - u) x w, 



with Tf, = a 2 / (6^) the characteristic bubble response time, and u and oj denoting re- 
spectively the fluid velocity and vorticity vectors at the bubble position. Fully-elastic 
boundary conditions are assumed for the bubbles on the cylinder's lateral walls, while 
direct bubble-bubble interactions (4-way coupling) are neglected. It is important to note 
that the lift force in (|2.8[) is only approximate. Here the expression for the rot ational 
flow in the limit Ret — > oo is used, corresponding to a lift coefficient Cl = 1/2 ( AutonI 
(1987)). iowever, in the simulation we consider the case of Reb ~ O(10) in an unsteady 
flow for which the Auton's ex pression may be i ncorre ct or may overestimate the ac- 
tual intensity of the lift force (jvan Nierop et al. (|2007h ). Despite the above mentioned 



ap proximations and the sim plifie d approach of the mode l adopted, it has been shown 
bv Ivan den Berg et al. ( 20061 ) and Calzavarini et al. ( 20061 ). that a qualitative agreement 



with real systems occurs. We will further discuss on the effect of the lift on the bubble 
dispersion in the TC flow in section l3~5l of this paper. 

2.3. Numerical implementation, code validation and tuning of the DNS 

The continuous phase, (|2.5[) . has been implemented in cylindrical coordinates on a finite- 
difference scheme. The grid, N z x N r x Nq = 128 x 64 x 256, is non uniform and refined in 
the near-wall regions in the radial direction, while it is uniform both in the vertical and 
azimuthal ones. Careful attention has been paid to the implementation of the non- linear 
term by the use of an algorithm that enjoys highly-conservative properties for energy, 
proposed bv iFukagata fc Kasa gi (2002). The time marching scheme is of second order 
accuracy. The bu bble trajectories (12 7I) - (|2~8| have also been implemented in cylindrical 
coordinates as in lDieridi et al\ ([2004) and integrated by a second order scheme in time, 
while interpolation of fluid velocity at the bubble position is implemented by a third- 
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order scheme, i.e. taking into account up to the second nearest lattice nodes (64 points). 
The feed-back ()2.6p is extrapolated from the bubble position to the next nearest nodes 
by a tri- linear scheme (8 points). We use periodic boundary conditions in the vertical 
direction both for the fluid and the bubbles. Additionally, for the fluid only, we fix the 
mean flow in the vertical direction to be zero, in order to avoid the emerging of a mean 
axial current that would prevent a realistic comparison with the experiment, where it is 
forbidden for geometrical reasons. 

The implementation of the single-phase flow dynamics has been validated by measuring 
the onset of the primary instability, i.e., the critical Reynolds number as a function of 
the gap ratio, Re c (r)), as well as the dimensionless torque G — T/(pv 2 L) at different 
Reynolds numbers and for different values of (the gap r atio) n . The results on Re c (77) 
were compared w ith the available analytic a l calc ulations Tavlor (|l923l ) (for the narrow- 
gap limi t 77 — > 1). iDominguez-Lerma et ~al\ (Il984h and with phen omenological estimates 
given bv lEsser fc Grossmannl (|l996f ) and lDutcher fc Mulierl ( 20071). Moreover we c ompa re 
G(rj,Re) with the empirical law given by IWendtT i 19331 ) "(also in Lathrop et al. ( 1992f )). 
finding satisfactory agreement (figure [2]). 

The tuning of our DNS with the experimental setup adopted by Murai et al. is then 
done as follows. The experimental apparatus is geometricall y characterized by a gap 
ratio 77 = 5/6 and an aspect ratio r = L/d = 20, see figure 1 of lMurai et all (|2005l ) for a 
sketch of their setup. The flow boundaries adopted in the experiment on the horizontal 
surfaces are free fluid surface (free-slip) at the top and zero mass-flux at the bottom. 
Since we do not expect drag to be strongly affected by the influence of the top/bottom 
boundary conditions, particularly for aspect ratio values (r) much larger than one, we 
use, for simplicity and in order to limit the computational costs a smaller aspect ratio 
Y/d = 8 than in experiment, and periodic boundary conditions as mentioned before. 
The gap ratio is as in the experiment, r\ = 5/6. It turned out that the constraint of 
periodic boundaries to the flow may introduce a bias on the formation and breaking of 
coherent structures in the flow, on this point we will comment more in detail in section 
3.2. Furthermore, for the whole set of numerical simulation we fix the ratio a/d and 
the bubble Reynolds number Reb to the experimental values - both these quantities are 
not Reynolds dependent in the experiment - and we vary the fluid Reynolds number 
Re, by changing the inner cylinder rotational speed Q as in the experimental setup, in 
the range Re = 600-2500. Care has to be put on the tuning of the void fraction, since 
in the experiment it was not kept constant among the different runs but estimated a 
posteriori through a measure of the total gas flow rate. We chose the total void fraction 
cto consistent with the experiment. The values (Tabl e [J) a re obtained by a fit to the 
estimated experimental void fraction (see also lOiwal (|2005l) . chapter 5); this explains 
some tiny deviations among the two cases. For the flow conditions that we consider for 
the comparison, i.e. Re ^ 2500, ao is always less than 1%. The value of the ratio a/d 
and the total void fraction value ao, then fix the number of bubble we must use in the 
numerical study, namely Nf, — ao(TrL(R 2 — Rf))/^ira 3 . The number of bubbles Nt, varies 
between few thousands to a maximum of roughly 1.5 • 10 4 . Table [T] summarizes the set 
of parameters adopted for the simulations and the experiment. 

A typical simulation run at a given value of the Reynolds number is composed of two 
parts. First, a single-phase DNS is performed until statistically stationarity is reached. 
Then bubbles are positioned in the flow with random homogeneous distribution and ve- 
locities equal to the local fluid velocity. This bubble injection procedure is supposed to 
mimic the presumably fully developed state at the top TC section where the bubble 
distributions have been measured in experiment. Simultaneously with the bubble injec- 
tion, the bubbles-fluid (two-way) coupling is activated. The two-phase simulation is then 
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| V | r | a/d 



b.c. 



Re 



I ao (%) 



\Re b \ We 



Exp. 



5/6 



20 



1/40 



vert, no-slip, 
horiz. top free-slip, 
horiz. bottom zero-mass flux 



600 - 4500 



0.092 at Re = 600 
0.721 at Re = 2700 
2.446 at Re = 4200 



7.1 



sC 0.6 



Sim. 



5/6 



1/40 



vert, no-slip, 
horiz. periodic 



600 - 2500 



0.125 at Re = 600 
0.670 at Re = 2500 



7.1 



Table 1. Releva nt scales and param eters for the two-phase Taylor-Couette system studied 
experimentally bv lMurai et al\ (| 20051 ) and numerically in this paper. We report first the three 
geometrical parameters: gap ratio n, aspect ratio F, and the ratio between the bubble radius 
and the gap width a/d, then the flow boundary condition, and the large scale flow Reynolds 
number range. Finally, bubble characteristics are listed: the global volume void fraction ao (in 
per cent), and the bubble's typical Reynolds and Weber numbers. 




2500 



Figure 3. Ratio a/rj between bubble radius and Kolmogorov scale (circles) and ratio 
a + — a^/r w / (pv 2 ) between bubble radius and wall length unit, both as a function of Re. 
Both the Kolmogorov length n and the wall-shear-stress t w are evaluated on the single-phase 
flows. The thickness of the viscous-sub layer (l v ) in a planar channel-flow in wall-units is roughly 
L, — 5. 



performed until a statistically steady state is well established. All relevant physical ob- 
servables in the single- and two-phase flow are computed by suitable time averages, of 
order O(10 2 ) revolution times (= 27rf2~ 1 ), in the statistically stationary regime. 

Finally, in order to verify the consistency of the model adopted in the DNS, we check 
if the smallest scale in the system is of the same order of the bubble typical size. In 
particular we compare the bubble radius a with the size of the Kolmogorov dissipative 
scale 77 = (i/ 3 /(e)) 1 / 4 and, alternatively, we look at the magnitude of the wall-units 
radial size a + = a^J t w / {pv 2 ) . While the ratio a/77 is relevant in the bulk flow the a + 
is of importance in near-wall regions. The results reported in figure [3] show that our 
assumption of point bubbles indeed is reasonable. 
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Figure 4. Torque reduction ratio compu t ed in the numerical simulations and the experimental 
measurements performed bv lMurai et al\ (|2005l ). In the inset TRR is divided by the mean void 
fraction qq. Note that the void fraction is not kept constant with Re, see table [1] 



I 



2.5 



2.0 



1.5 



1.0 



0.5 



0.0 



fi£>=600 




Re=900 




ffe=1200 




Re 2000 




Re 250(1 






// 




Jf 









0.0 



0.2 




0.N 



1.0 



0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 

(r-Ryd {r-R<yd 

Figure 5. Normalized mean local void fraction in the radial direction, a(r)/cto, at different 
Re values, in the numerics (left) and the experiment (right). The quantity a(r)/ao has been 
adopted from the data in figure 9 of lMurai etal\ (|2005f >. For better comparison the homogeneous 
distribution line a(r) /qq = 1 is also shown. 



3. Numerical results 



3.1. Overview 



In this section we present the comparison of numerical and experimental results. In 
figure [4] we show both the TRR coefficient obtained from the experiment and from the 
simulations. I n the inset TRR is n ormalized by the total void fraction (this ratio is called 
sensitivity bv iMurai et al. ( 2005h ). The quantitative agreement on this observable is 
satisfactory. First, we detect a robust (up to 20%) statistically steady drag reduction. 



Microbubbly drag reduction in Taylor- Couette flow 



9 



Second, drag reduction effect has a decreasing trend at increasing the Reynolds number. 
How can the two above mentioned observations be explained? Paragraphs 13.21 and 13.41 
of this section will be devoted to understand the physical mechanism leading to this 
behavior. 

Less satisfactory is the comparison with the experiment of the local mean void fraction 
(or mean bubble concentration). In figure[5]we show the mean radial bubble concentration 
profiles normalized by the total void fraction, both from the numerical simulations and 
from the experiment. This quantity is defined as, 

a{r)_Rl-R* 1 f*°+ T £ ^ S(r y r[l) (t)) ^ 



a 2N b T s ] ta A ~ J r 

where y r (i)(t) is the radial component of the position vector of the i-th bubble and T s 
indicates the duration of the time interval where the two-phase flow is in statistically 
steady conditions. The local bubble concentration in the simulations is less inhomoge- 
neous than in the experiment. The maximum deviation from equi-distribution is about 
100% of the homogeneous case and is peaked both at the inner and the outer wall. On 
the other hand, in the experiment, where radial mean void fraction has been evaluated 
by averaging over a series of snapshots taken in the upper part of the TC system, a much 
stronger accumulation peaked near the inner cylinder is found and it reaches much larger 
values (roughly 4 times larger than in numerics in the most extreme case Re = 600). 
Nevertheless, the same trend towards homogenization at increasing Reynolds number is 
observed in both cases. The possible origin of this discrepancy and its relevance on the 
effect of drag reduction will be further discussed in detail and will be the focus of section 
331 



3.2. Origin of the drag reduction 

The variation of the intensity of the torque in the TC system is directly connected to a 
change in the mean azimuthal velocity profile along the radial direction ue(r), see (|2.2|) - 
(|2.3[) . Figure [5] (top) shows the actual modification of the two-phase profile compared to 
the single-phase at two different Reynolds numbers, Re = 900 and Re = 2000. These Re 
values will be adopted as reference cases in the following discussion. A smoother slope 
is observed in figure [5] for the gradient at the wall of Tlq (r ) in the two-phase case. As 
consequence the torque is reduced. We also note that the torque Tt for the two phase 
system, when normalized by the torque 7} = pfl 2 RfL/ + rf)Re) in the (single-phase) 
laminar regime, satisfies the following relation (derived in the appendix): 

where u r ug{r) is the azimuthal-radial component of the mean Reynolds shear stress 
and fb,e{r) the mean azimuthal component of bubble forcing. Relation (|3.2[) reveals the 
connection between the Reynolds shear stress, the bubble forcing, and the torque. We 
have tested that in the two-phase case the mean bubble forcing ff, is mainly dominated 
by the vertical component, see figure [7] Therefore in the relation (|3.2j) the contribution 
coming from fb,e(r) is negligible compared to the Reynolds shear stress and the measure of 
u r ug(r) is a direct way to appreciate the variations in the mean flow features responsible 
for a torque change. As shown in Fig. [6] (bottom), this variation appears to be relevant 
in the whole bulk of the system and is approximately constant through it, both at large 
and small Reynolds numbers. 

We also have previously observed that in the relative low Reynolds range under con- 
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Figure 6. Top panel: mean azimuthal velocity profile ug(r-) at Re — 900 (left) and Re = 2000 
(right). The laminar profile is also reported. Bottom panel: the mean Reynolds shear stress 
component ugu r (r) at Re = 900 (left) and Re = 2000 (right). 
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Figure 7. Normalized globally averaged (space and time) relative components of the bubble 
forcing \(fb,g)\/(\{fb,g)\ + \{ft,r)\ + |(/ M )|) with (3 = z,r,9 vs. Re. 



sideration, coherent flow structures in the form of Taylor vortices are persistent in the 
flow. To reveal the effect of the interaction of the dispersed bubbly phase with the co- 
herent flow structures, we present a visualization of the flow. Figure [8] shows snapshots 
from two simulations respectively at Re = 900 and Re = 2000, both in the single-phase 
and in two-phase case. Vortical structures are identified by iso-surfaces of the second 
invariant of the velocity gradient tensor and are shown in a uniform color accordingly to 
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Figure 8. Flow visualization at Re — 900 (top) and Re = 2000 (bottom): Single-phase flow (left) 
and two- phase flow (right). The Taylor vortex structures are identified by iso-surfaces defined 
by the relation (diUj)(djUi) = —0.4, and uniformly colored (white and violet) accordingly to 
the sign of too on the surfaces. The wall shear stress, measured in p(Rift) 2 units, is reproduced 
in color on the inner cylinder surface. 



the sign of the azimuth al component of the local vorticity field, ue 



Tanahashi et al 



( Hunt et all (|1988f ). 

( 200ll )). The intensity of the inner-shear stress, T W i, is also reproduced 
in colors in the figure. Vortical Taylor structures, that come in counter-rotating pairs, are 
well captured by this visualization method. In both cases it is evident how the addiction 
of microbubbles efficiently breaks the self organized order established by the system in 
the single-phase condition. 



3.3. Effect of bubbles on Taylor vortices 

We now have a more detail look at the Taylor vortices and how they are affected by the 
bubble injection. We note that in a periodic TC flow Taylor vortex rolls may appear only 
in even numbers. A Fourier spectral analysis of the azimuthal vorticity component on 
the mid-cylindrical surface of the gap, uig(r = Ri + d/2, z), reveals that the numerically 
simulated single-phase TC flow displays respectively 8 Taylor vortices for Reynolds values 
Re ^ 2000 and 6 above this threshold (Re > 2000). The two-phase flow on the other 
hand shows a different behavior at small Reynolds number, in particular for Re ^ 900 
we observe only 6 Taylor vortices, while the same number of rolls as the single-phase is 
realized for Re > 900. This is qualitatively consistent with the experimental observation 
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made bv lMurai et al\ (|2005h . as they observe an elongation of Taylor vortex spacing (i.e. 



the center-to-center distance) in the two-phase case as compared to the single-phase case. 
Such elongation decreases as Re is increased. The single-phase/two-phase elongation ratio 
in the low-Reynolds number regime Re ^ 900 is roughly 30% both in the simulations and 
in the experiment. However, we cannot rule-out that the even-parity constrain imposed 
by the periodic boundary condition in the DNS may introduce some kind of artifact 
on the overall st ability of Tay l or vo rtex couples. Also, we note that the experimental 
measurement, see lMurai et al. I d2005h figure 13, has been performed with larger bubbles 



(a = 1mm), as compared to the TRR measurements discussed in that paper (a = 0.3mm). 
The bubbles employed to quantify the vortex elongation in the experiment turn out to 
corresp o nd ro ughly to i?e b = 400 and We = 90 and they are normally in a wobbling state 
(jCxracd (ll973h V Such conditions lie outside the range of applicability of the Eulerian- 



Lagrangian model and they prevent from attempting a one-to-one comparison with the 
DNS. 

We now focus on the azimuthal wave-length dependence in the DNS flow. We note 
that any signature of a dominant wave-length present in the single-phase case is removed 
from the spectra of uJe(r = R t + d/2, z) in the two-phase flow case. A transition from a 
clear azimuthal wave-length, that characterize the wavy-vortices and it is here of wave 
number 4, to a continuous and noisy-dominated spectra, is observed. 

Again from this analysis we realize that the bubble phase seems to strongly break 
the coherence of the self-organized vortical structures , as observed in th e previous flow 
visualizations (figured]). Therefore, differently from the Murai et al. (2005) interpretation 



given for the experiment, we associate the observed two-phase flow drag reduction in 
the DNS to the breaking and consequent re-arrangement of coherent flow structures 
rather than to its stretching or displacement. This view will be further supported by the 
discussion in the next subsection. 

3.4. Explanation for TRR decrease at increasing Re 

We now want to address the second of our questions, namely, why the drag reduction 
becomes less and less when increasing the inner cylinder rotational speed. Therefore we 
look at the global energy balance relation for the torque in the two-phase flow (see the 
appendix), 

T b = "* B Z-W L « £6 ) - <f b • u» (3.3) 

= «e b ) - <f b • u» . (3.4) 

The second equation (|3 .4|) follows from relation (|2.4p . TRR = 1 — T/T b is thus composed 
of the sum of the two terms 

which are here both of importance. While the former is associated to changes of the 
dissipative structures in the system, the latter is associated to the extra energy input in 
the system due to the bubble momentum transfer. Both these terms are positive (i.e., 
contribute to the drag reduction). Furthermore, the 1 — (e b )/(e) term is dominant at 
all Re, see figure |H This means that the bubbly phase, that breaks vortical dissipation 
regions, modifies the energy flux directed from the large-scale of the external forcing (due 
to the rotating cylinder) to the dissipative scale. 

Both contributions 1 — (eb)/(e) and (f b • u)/(e) to TRR are vanishing at increasing 



Microbubbly drag reduction in Taylor- Couette flow . 



13 



0.15 



0.10 



0.05 



0.00 




2500 



Figure 9. The correlation coefficient (fj, • u)/(e) and the relative variation of the mean energy 
dissipation rate of the two-phase flow as compared to the single-phase, 1 — (eb)/{e) as function 
of Re. The sum of the two terms is the total torque reduction ratio (TRR). 
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Figure 10. Mean vertical fluid velocity at the bubble positions, in dimensionless units 
(u z )b/(£lRi) in a one-way and two-way run at variable Re. Bubbles stay preferentially in 
down-flow regions in the former case while is just the opposite for the latter. 



Re. How to explain this dependence? We note that (ff, • u), the correlation term of 
bubble forcing with velocity fluctuations, can be positive for two different reasons. First, 
the bubbles might distribute preferentially in regions where the velocity is up-flow as 
compared to the mean (and for clarity we note that the global mean velocity in the TC 
system is zero). This would produce a positive correlation with the term — g that is the 
dominant one in the bubble feed-back f;,. The second possibility is that bubbles may 
induce a local up-flow fluctuation in the surrounding fluid that is larger as compared to 
the underlying flow fluctuation. This would produce always a positive correlation with the 
gravity term — g. The latter mechanism, where the bubble perturbation is overwhelming, 
is similar to what happens in a system where bubbles are rising in an otherwise quiescent 
fluid. Therefore, we will call this type of dynamics pseudo-turbulent. 

In other words, given the assumption (f;, • u) ~ g • (u)b, where (. . .) b denotes the 
average on the bubble centroids, we would like to understand the behavior of (u}(, in 
the different two-phase flow regimes. To distinguish between the two above mentioned 
possibilities, we make a numerical test by switching-off the two-way coupling, i. c., f(, = 0. 
Also in this one-way coupling case bubbles distribute inhomogeneously. In particular, we 
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Figure 11. The ratio A z /A z p versus Re, with A z = {u' z 2 )/{u' 2 +u' r 2 ) and u' = ((u - u(r)) 2 ) 1/2 
mean velocity fluctuation, evaluated for the two-phase flow Af p represents the same quantity 
evaluated in the single-phase case. In the inset, both A z and Af p are shown separately as a 
function of Re. 

expect that the bubble distribution in the two-way and one-way coupling cases shares 
similar properties if the amplitude of fluid fluctuations compared to the bubble feed- 
back is strong enough, while we expect them to have different properties if the pseudo- 
turbulence mechanism holds. The evaluation of the quantity (u z }bt m particular as a 
function of Re, allows to discriminate between the two cases. In figure [10] we show its 
behavior both in the one-way and two-way coupling case. In the one-way coupling case 
we observe preferentially bubble accumulation in down-flow regions. More precisely at 
low-i?e, from flow and bubble visualization, we may observe that bubbles, while rising, 
accumulate on stable regions in the flow. The phenomenology of the large-i?e regime 
is instead rather different, since the flow is time-dependent we do not observe bubble 
accumulation in stable positions but a rather homogeneous distribution that slightly 
favors downflow-regions. If we look now to the two-way coupled case we realize that 
bubbles strongly modify the flow by producing an upward flow in the position where 
they are. Therefore this test supports the hypothesis of a pseudo-turbulent mechanism 
acting in the two-phase flow. The pseudo-turbulent mechanism must become less effective 
at large Re numbers because of the increased level of fluid fluctuations. This is confirmed 
by the measure of the intensity of the vertical anisotropic energy content in the system, 
figure QT] Vertical anisotropy in the two-phase flow, namely A z = (u^)/{ug + u r 2 ) with 
u' = ((u — u(r)) 2 } 1 / 2 mean velocity fluctuation, when normalized by the anisotropy in 
the single-phase case, , monotonically decreases at increasing the Reynolds number. 

To summarize, the mainly vertical extra-forcing exerted by the bubbles on the fluid is 
responsible for local and anisotropic modifications of the fluid fluctuations. No large-scale 
flow mechanism induced by bubbles seems to be present or relevant in this process. At 
increased external forcing (larger Re) , the efficiency of the bubbles in perturbing the flow 
is lost. 
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Figure 12. Normalized radial mean void fraction at Re = 600. Here a comparison of the 
numerical r esults with and wit hout lift force are shown, together with the experimental profiles 
reported in lMurai~e £ al. (2005|). In the case without lift, nearly all the bubbles collapse into the 
first compartment of the histogram, giving a sharp bubble concentrations whose peak has been 
truncated for better readability. 



3.5. Bubble distribution 

Inhomogeneous bubble distribution in the TC system is due to the concurrence of several 
competing factors: (i) the centrifugal force, whose intensity towards the inner wall is 
proportional to —u„/r; (ii) inert i a tha t pushes the bubbles towards the local vortical 
flow structures (jWang fc Maxevl (jl993l )); (hi) th e drift towards the w alls, that is a lift 
effect due to the presence of a local mean shear (jSerizawa et all (jl975l ); and finally, (iv) 
the pseudo-turbul ence flow fluctuati o ns which may produce a wea k side by side bubble- 
bubble attraction (jZenit et al. ( 2001 ). Bunner fc Trvggvasonl ( 2002h ). All these effects are 
included explicitly, or via the bubble-flow interactions, in the numerical model system. 

For simplicity we limit our investigation to the radial dependence of mean bubble con- 
centration profiles. The formations of vertic al columnar or spira ling structures, similarly 
to the ones experimentally investigated bv lShiomi et al\ (|l993[ ), are indeed observed in 
our numerics, but this will not be the focus of this section. Already from numerical results 
reported in figure Owe have observed that at low Re numbers bubbles accumulation at 
both walls is observed, while at large Re a weak tendency of accumulation in the core 
is noticeable. Since in the low Re condition the drag behavior is dominated by pseudo- 
turbulence, we attribute near-wall accumulation to the bubble- wall interaction due to the 
lift in our numerical simulations. More explicitly, the local vertical mean flow produced 
by the bubble coupling is constrained by the no-slip conditions at t he walls. This local 



up-flow profile, as in an up- flow vertical pipe ( Serizawa et al. ( 19751 )). induces a lateral 



(lift-driven) migration to the walls. Indeed this is evident if we switch off the lift force. 
Without lift all the bubbles collapse on the inner cylinder due to the centripetal force, 
see figure [TSJ At large Re the pseudo-turbulent effect is less efficient and the bubbles 
tend to accumulate in vortex core regions: inertia plays the dominant role. Figure [T2l also 
shows a comparison at Re = 600 with the bubble concentration profile measured in the 
experiment. As noted previously, in the experimental profiles the local void fraction is 
always peaked near to the inner cylinder, in contrast to what is found in the simulations. 

What is the main reason for this discrepancy? First, the difference in statistically 
stationarity of the bubbly phase distribution. In the numerics the bubbles are injected 
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statistically homogeneously into the fluid domain and a relatively long time, typically 
0(100) revolution times, is required until statistically stationarity for the bubbles dis- 
tribution is reached. In particular, after the bubbles release in the flow, we observe a 
transient behavior in their distribution. Bubbles first accumulate in separate columnar 
regions and then, when flow structures are further broken, they merge in larger domains. 

Next, as we have already noted in the introduction, different bubble injection proce- 
dures are employed in the numerics as compared to the experiments. There, injection 
sites are localized at the bottom and near to the inner cylinder wall. From that positions 
bubbles spread to the upper section of the TC flow, where the bubble distributions are 
finally evaluated. So, it could be that under the particular conditions of slow rotation of 
the inner cylinder (low Re) stationarity is not fully achieved in the experimental set-up, 
since the bu bbles rise almost ver tically. This seems to be supported by some experimental 
snapshots in iMurai et all (|2005h . see figure 8 of that paper and the corresponding discus- 
sion. Furthermore, bubble-bubble coalescence may sporadically occur in the experiment. 
Coalesced bubbles, which rise faster than the others, lead to the formation of clusters 
and may delay or prevent the bubble distributions to become statistically stationary (Y. 
Murai 2007, personal communication). 

Second, on the numerical side, a possible source for the observed discrepancy may lie 
in the incompleteness of the model adopted for the lift force that may overestimate its 
relevance. It has been shown (jMagnaudet fc Legendrd (|l998T that for non-deformable 
bubbles, whose Reynolds number (Reb) is larger than one, and in a pure shear flow the lift 
coefficient can be considerably smaller than the asymptotic valu e, Cl = 1/2. Similar find - 
ings for the lift coefficient have experimentally been obtained by lvan Nierop et all (J2007J). 
Moreover, Cl depends on the local shear and the local vorticity. In order to understand 
how robust the shape of the mean bubble concentration profile is when changing the lift 
coefficient Cl, we have performed a test at the lowest Reynolds number (Re — 600), 
where the effect of bubble accumulation on the wall is more intense. Cl has been varied 
in the interval [0, 0.5] for a set of numerical runs all starting from the same initial con- 
dition and extended until statistically stationarity in time is reached. As shown in Fig. 
[T51 it turns out that the bubbly mean concentration is rather sensitive to Cl- For the 
cases with Cl ^ 0.1 it exhibits a single wall-peak near the inner wall. This also implies a 
dependence on the bubble feed-back on the fluid. As reported in figure [T3l (inset), reduc- 
ing the lift coefficient progressively from Cl = 0.5 to zero corresponds to a reduction in 
TRR. In particular, the accumulation of a small volume concentration of bubbles (here 
«o = 0.125%) on the surface of the inner cylinder, as for the cases Cl — and Cl — 0.1, 
would have almost no effect on the torque (or drag) modulation. 



4. Conclusions 

We have numerically investigated the effect of microbubbles on a Taylor-Couette (TC) 
flow in the wavy vortex flow regime and focused on a comparison with a specific ex- 
periment where variations in the drag have been observed. We demonstrated that this 
phenomenon is reproduced by an Eulerian-Lagrangian model with point-force coupling. 
The main effect of microbubbles is to create a local perturbation of the flow, mainly di- 
rected upwards along the vertical direction that is able to break the coherent, and mainly 
dissipative, vortical structures of the flow. This pseudo-turbulent mechanism is at the 
origin of the torque reduction in the TC flow regime considered here. As a consequence at 
large Reynolds number (Re), where the typical fluid velocity fluctuations are larger and 
coherent structures are unsteady and less persistent, this mechanism loses its efficiency. 
The lift force resulted to be particularly important in this process. The strength of the 
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Figure 13. Effect of the intensity of the lift force on the bubble mean radial distribution and on 
the drag reduction. The lift coefficient Cl is varied in the interval [0, 0.5]. The left panel shows 
a(r)/oto at changing Cl- On the right the corresponding shape of the mean azimuthal velocity 
profile, ug(r), is shown. The inset reports the corresponding value of TRR versus Cl. 



lift force is the main factor for the bubble mean concentration profiles. Furthermore, the 
drag reduction is particularly sensitive to it. When suppressing the lift, almost no drag 
reduction is obtained. 

The dynamics responsible for the drag reduction in the low Re TC flow highlighted 
in this study is fundamentally different from the one that has been observed in p revious 
highly turbulent experimental investigations by Ivan den Berg et al. (|2005l l2007h . From 
our numerics it can be extrapolated that small non-deformable bubbles (zero Weber 
number (We)), whose size is always below the smallest typical scale of the flow will be 
of less and less importance for drag reduction as the external forcing (and henc e Re) is 
increased. This is also consistent with the theory developed bv lL'vov et al. ( 2005t ). where 
for small bubble in a channel flow the variation of the drag expected i s only of the order of 



the void frac tion. This suggests that at larger Re, when We > 1 (as in Ivan den Berg et 



2005L 2007)), o ther physical mec hanisms shall b e effective, na mely, bubble deformatio n 
Lu et a/.l (|2005l )) or compression (|Lo et al\ (|2006h ) or splitting (|Meng fc Uhlmanl (|l998h ). 



Further joint numerical and experimental works are still needed to clarify this aspect. 
Acknowledgment: We acknowledge Dr. Y. Murai for giving us some essential informations 
concerning his experimental setup and measurements and for reading the manuscript. 



Appendix A. 

In this appendix we present a derivation of the equation s (|3.2[) and (|3.3p . based on the 
use of the reciprocal theorem ( Happel &: Brenner! (1973)). Let us introduce two incom- 
pressible velocity fields u and u, both satisfying the same set of boundary conditions. 
We identify u as the real flow in the TC system, and u as a particular flow field, that 
will be specified afterwards. We focus here on the axially periodic TC system with only 
rotating inner cylinder. Therefore, the wall boundary conditions are the following: 



u = u = Riflee at r = Ri 
u = u = at r = R 



(Al) 
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The instantaneous, i.e. time dependent, torque T acting on the inner cylinder is given 

by 



(A 2) 



TQ, = R i fl(f> d 2 x (-e r ■ er ■ eg), 

Jr=Ri 

where cr is the stress tensor associated to u, and e r and eg are unit vectors in the radial 
and azimuthal directions, respectively. Using the boundary conditions (|A ip for u, we can 
rewrite (|A 2[) in a volume integral form 



Ttt= <L d 2 x n • a ■ u = / d 3 x (V • a) ■ u + / d 3 x er : Vu. 

/s Jv Jv 



(A3) 



Here £> is the area of the boundaries, V is the fluid v olume bounded by S , and t he double- 
dot represents a product of two dyadics (see e.g. lHappel fc Brenner! (|1973T ). Sec. 2-1). 
Substituting the momentum equation 



V-cr = p(a t U+(u- V)u-f 6 ), 



into (IA 3D. we obtain 



TQ = p d 3 x (d t u) ■ u - (uu) : S f b ■ u 



d 2 x (n • u)(u • u) + / d 3 x er : Vu, 

=o 

(A 4) 

where S(= {Vu+ (Vu) T }/2) denotes a strain tensor. Similarly to the derivation of (|A 3j) . 
one can obtain the inner torque T in the field (u, cr) 



(A 5) 



Tfl = / d J x (V • or) • u + / d 3 x cr : Vu 



We no w apply the Lorentz reciprocal theorem cr : Vu = e> : Vu (see e.g. lHappel fc Brenner 
(Il973l) Sec. 3-5) to ([X4| and jX5j, and rewrite the difference of the two torques in the 
volume integral form, 



(T - t)fi = p / d 3 x (9 t u) ■ u - (uu) : S - f b • u 



d 3 x (V • or) ■ u. 



(A 6) 



From equation (|A 4|) and (|A 6p the relevant relations (|3.2p and (|3.3p can readily be ob- 
tained by proper assumptions on the field u. We consider two cases: (i) u is the real flow 
field u, (ii) u is the steady flow field solution of the Stokes equation for the TC system. 

Case (i): If we choose the field (u, <x) to be (u, cr), the product (cr : Vu) is equal to the 
energy dissipation rate e(= 2vS : S) and (|A 4j) becomes 



dt 



d 3 x 1 ' 



d x (n • u) 



pu ■ u 



d 3 x (e - f 6 • u) , 



(A 7) 



=o 



which corresponds to the kinetic energy transport budget. In a statistically steady state 
((^ . . .) = 0), we obtain 



np(Rl - Rj)L 



= <e) - (fb ■ u) 



(A! 



which accounts for the contribution of the energy dissipation rate and bubble forcing to 
the time averaged torque Xh = (T) . 
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Case (ii): We now assume that the field (u, or) obeys the Stokes equation 

V • or = 0. 



(A 9) 



Imp osing the boundary co nditions (|A 1[) . we write the solution of circular Couette flow 
(see IChandrasekharj (|1961l) . Chapter VII) 



eo 



-Rix 



ri - m 



-(e r ee + e e e r )^^(RiQ). (A 10) 



Substituting (|A~9)) and (|A 10|1 into 
the mean torque 



Rl-Ri 

61) , and assuming a statistically state, we determine 



n = t, i 



Re{R 2 Q - Rj) V I u r u e (R 2 - r 2 )f b j 



2Rh 



(All) 



(A 12) 



2(i?^) 2 (l-77) 

where T; denotes the torque in the laminar flow, 

_ ^pjR^fRjL 
1 i?e 77(1 + 77) 

The first term inside the brackets on the right hand side of (|A lip accounts for the 
laminar flow contribution and the second term for the modification from the circular 
Couette flow due to the nonlinearity of the fluid motion and the bubble's forcing. Using 
the area averaged quantities u r ue(r) and fb,e(f), we may finally write (|A lip in the radial 
integral form, previously used in this paper, 

i-Ra 



T b = Ti\l 



Re 77 



(i?^) 2 (l-T7) 



dr 



Ri 



u r u e {r) _ (R 2 - r 2 )f bfi {r) ' 
r 2i?2 



(A 13) 



We note that a general identity linking the effect o f the Reynolds shear st ress distribution 
on the skin friction has been recently derived bv lFukagata et al. ( 2002h for the case of 
a pressure driven channel and circular pipe flows. The Fukagata et al. identity has been 
further extended to the probl em of channel flow with mixe d boundary conditions by 
using the reciprocal theorem in Sbragaglia fe; Sugivama ( 20071 ). Our derivation presented 
above represents an extension of that work to the TC system. 
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